#!/usr/bin/perl -w
use strict;

use FindBin;
use lib ("$FindBin::Bin/..", "/net/cpp-group/Leo/bin", "$FindBin::Bin");
use db_parameters;
use ortho_species;
use ensembl_parameters;
use chromosomes;
use min_max;
use die_error;
use run_cmd;
use DBI;
use align_gaps;
use Data::Dump qw(dump);

print STDERR <<"HEADLINE";
888888888888888888888888888888888888888888888888888888888888888888888

    table_coding_exons
	

	Description:
		Obtain the coding intron/exon structure of each transcript
	Populates DB tables:
		exon_loci
		prot_introns
		prot_exons
		prot_exon_stats
		prot_coding_exons

888888888888888888888888888888888888888888888888888888888888888888888
HEADLINE

use constant GENEID =>0;
use constant CHRMSM =>1;
use constant START_ =>2;
use constant END___ =>3;
use constant STRAND =>4;
use constant CENTRE =>5;
use constant ORDER_ =>5;


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

#	Determining exon structures of coding sequence


#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


#_________________________________________________________________________________________

#	retrieve_exon_coordinates

#_________________________________________________________________________________________
sub retrieve_exon_coordinates($$$)
{
	my ($dbh, $taxon, $chrm) = @_;
	print STDERR "Exon coordinates\n";
	return $dbh->selectall_arrayref(<<"PL/SQLCMD");
	SELECT
        DISTINCT ON (prot_id, prot_exons.start)
		prot_id,
		prot_exons.start,
		prot_exons.finish,
		strand,
        exon_id
	FROM
		taxon.prot_exons NATURAL JOIN
        ens_id JOIN 
		taxon.gene_loci USING (gene_id)
	WHERE
		chromosome = '$chrm' and
        species = '$taxon'
	ORDER BY 
		prot_id, 
		prot_exons.start;
PL/SQLCMD
}	

#_________________________________________________________________________________________

#	retrieve_cdna_sequences

#_________________________________________________________________________________________
sub retrieve_exonic_sequences($$$\%)
{
	my ($dbh, $taxon, $chrm, $cdna) = @_;
	my $cdna_data = $dbh->selectall_arrayref(<<"PL/SQLCMD");
	SELECT
		prot_id,
		exonic_dna_sequence
	FROM
        taxon.gene_loci NATURAL JOIN 
		ens_id NATURAL JOIN
		taxon.prot_exon_sequences
	WHERE
		chromosome = '$chrm' and
        species = '$taxon';
PL/SQLCMD
	for my $line(@$cdna_data)
	{
		my ($prot_id,  $cdna_seq) = @$line;
		$cdna->{"$prot_id"} = $cdna_seq;
	}

}

#_________________________________________________________________________________________

#	retrieve_cdna_sequences

#_________________________________________________________________________________________
sub retrieve_cdna_sequences($$$\%)
{
	my ($dbh, $taxon, $chrm, $cdna) = @_;
	my $cdna_data = $dbh->selectall_arrayref(<<"PL/SQLCMD");
	SELECT
		prot_id,
		cdna
	FROM
        taxon.gene_loci NATURAL JOIN 
		ens_id NATURAL JOIN
		cdna_sequences
	WHERE
		chromosome = '$chrm' and
        species = '$taxon';
PL/SQLCMD
	for my $line(@$cdna_data)
	{
		my ($prot_id,  $cdna_seq) = @$line;
		$cdna->{"$prot_id"} = $cdna_seq;
	}

}

#_________________________________________________________________________________________

#	map_exons_for_prot

#		map coding sequence onto the exons for this protein

#		also copies prot_coding exons to panda

#_________________________________________________________________________________________
use constant EXON_ID	=> 0;
use constant EXON_LEN	=> 1;
use constant EXON_LOCUS	=> 2;
use constant EXON_STRAND=> 3;
use constant EXON_START => 4;
use constant EXON_FINISH=> 5;
my $FH_exonic_mismatch_errors;
sub map_exons_for_prot($$\%\@\@\%)
{
	my ($dbh, $prot_id, $mapping_gaps_per_prot_id, $exons, 
		$exon_coding_structure, $prot_id_to_coding_exon_loci) = @_;

	#print STDERR dump($mapping_gaps_per_prot_id), "\n$prot_id\n";
	#ignore those which fail to map
	return unless (exists $mapping_gaps_per_prot_id->{$prot_id});

	# returns length of cDNA mapped onto exonal sequences.
	# i.e. $len may be longer than the length cDNA
	my ($beg_pos, $len, $gaps_str1, $gaps_str2, @gaps) = @{$mapping_gaps_per_prot_id->{$prot_id}};
	my $end_pos = $beg_pos + $len;

	#
	#	Failed: should print diagnositic to error file
	#
	unless ($len > 0)
	{
		print $FH_exonic_mismatch_errors
			  "$prot_id\n";
		return;
	}

	# holds exon_id, offset and coding length: 
	#		offset is 0 for all but first exon
	# 		length = exon length for all but last exon
	my @coding_struct;
	my @coding_struct_str;



	my $exon_cnt = 0;
	my $curr_pos = 0;
	my $exon_coding_start = 0;

	#print STDERR dump @$exons;
	#print STDERR join (", ",  $beg_pos, $len, $gaps_str1, $gaps_str2), "\n";
	
	#
	#
	#	Find first exon where coding begins
	#
	while ($exon_cnt < scalar @$exons)
	{

		#print STDERR "curr_pos = $curr_pos, beg_pos = $beg_pos, exon_len = $exons->[$exon_cnt][EXON_LEN]\n";
		#
		#	keep going if not within first used exon
		#
		if ($curr_pos + $exons->[$exon_cnt][EXON_LEN] <= $beg_pos)
		{
			$curr_pos += $exons->[$exon_cnt][EXON_LEN];
			$exon_cnt++;
		}

		#
		#	First used exon: Remember offset to exon start
		#
		else
		{
			$exon_coding_start = $beg_pos - $curr_pos;
			last;
		}
	}

	#
	# beginning not found in any exon
	#
	die "$prot_id\n" if ($exon_cnt == scalar @$exons);


	#
	#	
	#
	my ($genome_locus_beg,  $genome_locus_end);
	my (@exon_boundaries);
	#
	#	Find last exon of coding
	#
	while ($exon_cnt < scalar @$exons)
	{
		my $this_exon = $exons->[$exon_cnt];
		my $strand	= $this_exon->[EXON_STRAND];
		my $locus	= $this_exon->[EXON_LOCUS];


		#
		#	not last exon (Could be first exon)
		#
		if ($curr_pos + $this_exon->[EXON_LEN] < $end_pos)
		{
			my $coding_len_of_exon = $this_exon->[EXON_LEN] - $exon_coding_start;

			# save the genomic span for all exons
			$genome_locus_beg = $locus + $exon_coding_start * $strand unless ($genome_locus_beg);
			$genome_locus_end = $locus + $coding_len_of_exon * $strand;

			# save exon, offset and coding length
			push(@coding_struct_str, "$this_exon->[EXON_ID]:$exon_coding_start,". 
									$coding_len_of_exon);
			push(@coding_struct, [	$this_exon->[EXON_ID], 
									$exon_coding_start,
									$coding_len_of_exon]);
			
			# save exon loci
			push(@{$prot_id_to_coding_exon_loci->{$prot_id}}, 
				[$this_exon->[EXON_START], $this_exon->[EXON_FINISH]]);


			# increment by exon length
			$curr_pos += $this_exon->[EXON_LEN];
			$exon_cnt++;

			# save the exon boundary: 
			# i.e. the position of the end of the exon in the concatenated exonic DNA 
			push(@exon_boundaries,  $curr_pos);

			# only the first exon has an offset
			$exon_coding_start = 0;
			next;
		}

		#
		#	last exon (Could be only exon)
		#
		my $last_len = $end_pos - $curr_pos - $exon_coding_start;
		if ($last_len > 0)
		{
			# save the genomic span for all exons
			$genome_locus_beg = $locus + $exon_coding_start * $strand unless ($genome_locus_beg);
			$genome_locus_end = $locus + $last_len * $strand;

			# save exon, offset and coding length
			push(@coding_struct_str, "$this_exon->[EXON_ID]:$exon_coding_start,$last_len");
			push(@coding_struct, [	$this_exon->[EXON_ID], 
									$exon_coding_start,
									$last_len]);
			
			# save exon loci
			push(@{$prot_id_to_coding_exon_loci->{$prot_id}}, 
				[$this_exon->[EXON_START], $this_exon->[EXON_FINISH]]);

			# increment by exon length
			$curr_pos += $last_len + $exon_coding_start;
		}
		last;
	}


	#
	#	sanity checks
	#
	die "$prot_id,$curr_pos < $end_pos\n".dump(@coding_struct)."\n".dump(@$exons).
	dump(@{$mapping_gaps_per_prot_id->{$prot_id}}) if ($curr_pos < $end_pos);
	die unless defined $genome_locus_beg;
	die unless defined $genome_locus_end;

	# make sure genome coordinates are always beg < end
	if ($genome_locus_end < $genome_locus_beg)
	{
		($genome_locus_beg, $genome_locus_end) = ($genome_locus_end, $genome_locus_beg);
	}

	# get exon boundaries mapped onto cDNA
	my $exon_boundaries =  join (",", map { map_pos_via_aligned($_, @{$gaps[1]}, @{$gaps[0]}) } 
									@exon_boundaries);
#	my $exon_boundaries =  join (",", @exon_boundaries);
#	
	push(@$exon_coding_structure, join ("\t", 
                                        $prot_id,
										join ("|", @coding_struct_str),
										$genome_locus_beg,
										$genome_locus_end,
										$gaps_str1,
										$gaps_str2,
										$exon_boundaries)."\n");

	$dbh->do("COPY taxon.prot_coding_exons".
				"(prot_id, exon_id, coding_offset, length, ".
				 "rank, frame) FROM STDIN");
	my $cummulative_length = 0;
	my $frame = 0;
	for my $rank(0..$#coding_struct)
	{
		my ($exon_id, $coding_offset, $length) = @{$coding_struct[$rank]};
		$dbh->pg_putline(join ("\t", $prot_id, @{$coding_struct[$rank]}, $rank, $frame)."\n");
		$frame = ($coding_struct[$rank][2]) % 3;
		$cummulative_length += $length;
	}
	$dbh->pg_putline("\\.\n");
	$dbh->pg_endcopy;

}

#_________________________________________________________________________________________

#	map_cdna_onto_genomic_sequence_per_chromosome

#		for each transcript, provides mapping from exonic dna to cdna

#_________________________________________________________________________________________

sub map_cdna_onto_genomic_sequence_per_chromosome($$$\%)
{
	my ($dbh, $taxon, $chrm, $mapping_gaps_per_prot_id) = @_;

	#
	#	get sequence
	#
	my (%cdna,  %exonic_dna);
	retrieve_cdna_sequences($dbh, $taxon, $chrm, %cdna);
	retrieve_exonic_sequences($dbh, $taxon, $chrm, %exonic_dna);

	print STDERR "\t\tMapping chromosome $chrm: (", 
					scalar keys %exonic_dna, " transcripts: ";

	open_or_die (*TMPFA, ">$dir_pipeline_temp/delete.this.fa");
	
	my %transcripts_align_dp;
	my $cnt_exact_matches = 0;
	for my $prot_id (keys %exonic_dna)
	{
		die $prot_id unless exists $cdna{$prot_id};
		my $cdna_seq = $cdna{$prot_id};
		my $exon_seq = $exonic_dna{$prot_id};
		my $cdna_len = length($cdna_seq);
		#
		#	check if whole length match: do not need to do dynamic programming
		#
		if ($exon_seq =~ /^(.*)$cdna_seq(.*)$/)
		{
			my $initial_gap			= length($1);
			my $final_gap_offset	= $initial_gap + $cdna_len;
			my $final_gap			= length($2);
			my @gaps;
			# initial gap
			push(@gaps,  [0, $initial_gap]) if ($initial_gap);
			# end gap
			push(@gaps,  [$final_gap_offset, $final_gap]) if ($final_gap);
			$mapping_gaps_per_prot_id->{$prot_id} = [length($1), $cdna_len, 
										gaps_to_gapstr(@gaps),
										'',
										[@gaps],
										[]];
			++$cnt_exact_matches;
		}
		else
		{
			print TMPFA join ("\n", 
									'>'.$prot_id, 
										$cdna_seq, 
									'>'.$prot_id.'_exons', 
										$exon_seq, '');
			++$transcripts_align_dp{$prot_id};
		}
	}
	close TMPFA;
	print STDERR $cnt_exact_matches, " exact matches / ";
	my $cnt_dp_matches = scalar keys %transcripts_align_dp;
	print STDERR $cnt_dp_matches, " for gapped mapping / ";

	#print STDERR "\n";
	my $cmd = <<"CMD";
	seq_pairs_align -m dna < $dir_pipeline_temp/delete.this.fa
CMD
	if ($cnt_dp_matches)
	{
		my @results = get_results_run_cmd($cmd);
		for my $line (@results)
		{
			chomp($line);
			my @data = split /\t/, $line, -1;
			my $prot_id = $data[0];
			my $len_cdna = $data[7];
			my @gaps = (gapstr_to_gaps($data[-2]), 
						gapstr_to_gaps($data[-1]));
			my $beg2 = map_pos_via_aligned(0, @{$gaps[0]}, @{$gaps[1]});
			my $len2 = map_pos_via_aligned($len_cdna - 1, @{$gaps[0]}, @{$gaps[1]}) + 1 - $beg2;
			$mapping_gaps_per_prot_id->{$prot_id} = [$beg2, $len2, $data[-2], $data[-1], @gaps];
			delete $transcripts_align_dp{$prot_id};
		}
	}

	my $cnt_failed = keys %transcripts_align_dp;
	if ($cnt_failed)
	{
		print STDERR " $cnt_failed failed to map / " ;
		for my $prot_id(keys %transcripts_align_dp)
		{
			print $FH_exonic_mismatch_errors $prot_id, "\n";
		}
		#exit(1);
	}
	unlink("$dir_pipeline_temp/delete.this.fa");
	
}


sub save_prot_coding_structure($$$\%)
{
	my ($dbh, $taxon, $taxon_data, $prot_id_to_coding_exon_loci) = @_;
	print STDERR "\t$taxon Exons\n";

	#
	#	iterate chromosome at a time
	#
	my @chromosomes = get_chromosomes_with_sequences($dbh, $taxon, "$taxon_data->[SEQ_DIR]/dna");
#DEBUGG
	#@chromosomes = (2,3,4,5,6,7,8);
	print STDERR "\t\t", scalar @chromosomes, " chromosomes for $taxon\n";

	# save which exons contribute to the coding
	my @exon_coding_struct;
	
	for my $chrm(@chromosomes)
	{
		#
		#	mapping gaps
		#
		my %mapping_gaps_per_prot_id;
		map_cdna_onto_genomic_sequence_per_chromosome($dbh, $taxon, $chrm, %mapping_gaps_per_prot_id);

		#
		#	exon coordinates ordered by protein
		#
		my $exon_coordinates = retrieve_exon_coordinates($dbh, $taxon, $chrm);

		print STDERR "\t\tDeciphering coding structures (", 
						(scalar @$exon_coordinates), " exons / ";
		# for last exon of last protein
		push(@$exon_coordinates, ["LASTLAST"]) if @$exon_coordinates;


		my $cnt_seq = 0;
		my $curr_prot_id = '';

		# exon id and length making up current protein
		my @exons;
		
		# for each exon
		for my $line(@$exon_coordinates)
		{
			my ($new_prot_id,  $start, $finish, $strand, $exon_id) = @$line;

			if ($curr_prot_id ne $new_prot_id || $new_prot_id eq 'LASTLAST')
			{
				#
				# Start of new protein. Analyse exon structure of old protein
				#
				if ($curr_prot_id)
				{
					$cnt_seq++;
					map_exons_for_prot( $dbh, 
										$curr_prot_id, 
										%mapping_gaps_per_prot_id, 
										@exons, 
										@exon_coding_struct, 
										%$prot_id_to_coding_exon_loci);
					@exons = ();
				}
				last if ($new_prot_id eq 'LASTLAST');
				$curr_prot_id = $new_prot_id;
			}

			# get length of old protein
			my $len = $finish - $start + 1;

			#
			#	3' goes the other way
			#
			if (!$strand)
			{
				unshift(@exons, [$exon_id, $len, $finish, -1, $start, $finish]);
			}
			#
			#	Save 5'
			#
			else
			{
				push(@exons, [$exon_id, $len, $start, 1, $start, $finish]);
			}
		}
		print STDERR "$cnt_seq transcripts)\n";
		print STDERR "\t\tSaving Exons\n";

		$dbh->do("COPY taxon.prot_coding_structures".
					"(prot_id, coding_struct, start, finish, ".
					 "mapping_gaps1, mapping_gaps2, exon_boundaries) FROM STDIN");
		$dbh->pg_putline($_) for (@exon_coding_struct);
		$dbh->pg_putline("\\.\n");
		$dbh->pg_endcopy;
		@exon_coding_struct = ();

		#close $chromosome_FH;
	}
}








#_________________________________________________________________________________________

#	save_prot_coding_introns

#		calculates introns for all proteins in a species
#		also saves intron/exon counts per protein

#_________________________________________________________________________________________
use constant EXON_LOCUS_BEG =>0;
use constant EXON_LOCUS_END =>1;
sub save_prot_coding_introns($\%\@)
{
	my ($dbh, $prot_id_to_coding_exon_loci, $prot_exon_stats) = @_;

	print STDERR "\t\tSave Introns";
	$dbh->do("COPY taxon.prot_coding_introns".
				"(prot_id, rank, start, finish, length) FROM STDIN");
	while(my ($prot_id, $exon_loci) = each %$prot_id_to_coding_exon_loci)
	{
		#
		#	Get introns sorted by loci
		#
		@$exon_loci =  sort{$a->[EXON_LOCUS_BEG] <=> $b->[EXON_LOCUS_BEG] }(@$exon_loci);

		my ($intron_min, $intron_max) = (0,0);
		my ($cnt_exons, $cnt_disruptions) = (1,0);
		if (@$exon_loci >= 2)
		{
			($intron_min, $intron_max) = (1e99,0);
			for my $i (0..($#${exon_loci} - 1))
			{
				my ($start, 
					$finish) = ($exon_loci->[$i  ][EXON_LOCUS_END] + 1,
								$exon_loci->[$i+1][EXON_LOCUS_BEG] - 1);
				
				$dbh->func("$prot_id\t$i\t$start\t$finish\t".($finish - $start + 1)."\n",
							'putline');

				# save intron sizes for prot exon stats
				my $intron_len = $exon_loci->[$i+1][EXON_LOCUS_BEG] - $exon_loci->[$i][EXON_LOCUS_END] - 1;
				if ($intron_len > 10)
				{
					++$cnt_exons;
				}
				else
				{
					++$cnt_disruptions;
				}
				$intron_min = $intron_len if ($intron_len < $intron_min);
				$intron_max = $intron_len if ($intron_len > $intron_max);
			}
		}

		# save prot exon stats
		my $transcript_len = $exon_loci->[-1][EXON_LOCUS_END] - $exon_loci->[0][EXON_LOCUS_BEG] + 1;
		my $exon_count = @$exon_loci;
		push(@$prot_exon_stats, join ("\t", 
									  $prot_id, 
									  $intron_min, 
									  $intron_max, 
									  $exon_count, 
									  $transcript_len, 
									  $cnt_exons, 
									  $cnt_disruptions));
	}
	$dbh->pg_putline("\\.\n");
	$dbh->pg_endcopy;
}



#_________________________________________________________________________________________

#	save_prot_exon_stats

#		
#		saves intron exon counts, intron sizes, disruptions per protein

#_________________________________________________________________________________________
sub save_prot_exon_stats($\@)
{
	my ($dbh, $prot_exon_stats) = @_;
	#
	#	taxon.prot_exon_stats
	#
	print STDERR "/ Exon structure data to panda...\n";
	$dbh->do(<<"PL/SQLCMD");
    COPY
            taxon.prot_exon_stats
            (prot_id, 
             intron_min_size, 
             intron_max_size, 
             exon_count, 
             transcript_len, 
             real_exon_count, 
             disruptions_count) 
        FROM STDIN;
PL/SQLCMD
	$dbh->pg_putline($_."\n")  for (@$prot_exon_stats);
	$dbh->pg_putline("\\.\n");
	$dbh->pg_endcopy;

}



#_________________________________________________________________________________________

#	merge_disruptions_to_prot_exon_stats

#		
#		Merge disruption data into prot_exon_stats

#_________________________________________________________________________________________
sub merge_disruptions_to_prot_exon_stats($)
{
	my ($dbh) = @_;
	print STDERR "\tMerge disruption data into prot_exon_stats...\n";
	$dbh->do(<<"PL/SQLCMD");
    UPDATE
            taxon.prot_exon_stats
        SET
             disruptions_count = prot_disruptions.disruptions_count
        FROM 
            taxon.prot_disruptions
        WHERE
            prot_disruptions.prot_id = prot_exon_stats.prot_id;
PL/SQLCMD
}




#_________________________________________________________________________________________

#	save_coding_exon_counts_per_gene

#		
#		saves the number of coding exons per gene

#_________________________________________________________________________________________
sub save_coding_exon_counts_per_gene($)
{
	my ($dbh) = @_;
	print STDERR "\t\tSaving exon counts per gene to panda...\n";
	$dbh->do(<<"PL/SQLCMD");
	INSERT INTO
		taxon.gene_coding_exon_counts (gene_id, exon_count)
			SELECT
				gene_id,
				count(DISTINCT exon_id)
			FROM
				taxon.prot_coding_exons natural join
				ens_id
			GROUP BY gene_id;
PL/SQLCMD
	$dbh->do("ANALYSE taxon.gene_coding_exon_counts;");
}















#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

# Main logic

#88888888888888888888888888888888888888888888888888888888888888888888888888888888888888888



my $dbh = connect_to_panda();
use get_psql_time;
my $start_time = get_psql_time($dbh);

#
#	Get DB parameters
#
my %taxon_data = get_taxon_data;



#
#	taxon.prot_coding_exons
#

$dbh->{RaiseError} = 0;
$dbh->do("TRUNCATE taxon.prot_coding_structures") or die $dbh->errstr.".";
$dbh->do("VACUUM FULL taxon.prot_coding_structures;") or die $dbh->errstr.".";
$dbh->do("TRUNCATE taxon.prot_coding_exons") or die $dbh->errstr.".";
$dbh->do("VACUUM FULL taxon.prot_coding_exons;") or die $dbh->errstr.".";
$dbh->do("TRUNCATE taxon.prot_coding_introns;") or die $dbh->errstr.".";
$dbh->do("VACUUM FULL taxon.prot_coding_introns;") or die $dbh->errstr.".";
$dbh->do("TRUNCATE taxon.prot_exon_stats") or die $dbh->errstr.".";
$dbh->do("VACUUM FULL taxon.prot_exon_stats;") or die $dbh->errstr.".";
$dbh->do("TRUNCATE taxon.gene_coding_exon_counts;") or die $dbh->errstr.".";
$dbh->do("VACUUM FULL taxon.gene_coding_exon_counts;") or die $dbh->errstr.".";

#$dbh->do("DELETE FROM taxon.prot_coding_structures USING ens_id WHERE ens_id.prot_id = prot_coding_structures.prot_id and species = 'mono_ens'") or die $dbh->errstr.".";
#$dbh->do("VACUUM FULL taxon.prot_coding_structures;");
#$dbh->do("DELETE FROM taxon.prot_coding_exons USING ens_id WHERE ens_id.prot_id = prot_coding_exons.prot_id and species = 'mono_ens'") or die $dbh->errstr.".";
#$dbh->do("VACUUM FULL taxon.prot_coding_exons;");
#$dbh->do("DELETE FROM taxon.prot_coding_introns USING ens_id WHERE ens_id.prot_id = prot_coding_introns.prot_id and species = 'mono_ens'") or die $dbh->errstr.".";
#$dbh->do("VACUUM FULL taxon.prot_coding_introns;");
#$dbh->do("DELETE FROM taxon.prot_exon_stats USING ens_id WHERE ens_id.prot_id = prot_exon_stats.prot_id and species = 'mono_ens'") or die $dbh->errstr.".";
#$dbh->do("VACUUM FULL taxon.prot_exon_stats;");
#$dbh->do("DELETE FROM taxon.gene_coding_exon_counts USING ens_id WHERE ens_id.gene_id = gene_coding_exon_counts.gene_id and species = 'mono_ens'") or die $dbh->errstr.".";
#$dbh->do("VACUUM FULL taxon.gene_coding_exon_counts;");







for my $taxon(keys %taxon_data)
{
	my (%prot_id_to_coding_exon_loci);

	
	open_or_die(*FH_exonic_mismatch_errors, 
				">$dir_pipeline_errors/table_prot_coding_structures.$taxon.exon_genomic_mismatch.errors");
	$FH_exonic_mismatch_errors = *FH_exonic_mismatch_errors;

	#
	# Exons involved in coding per acc code
	#
	save_prot_coding_structure(	$dbh, 
								$taxon, 
								$taxon_data{$taxon}, 
								%prot_id_to_coding_exon_loci);


	my (@prot_exon_stats);
	#
	#	Introns involved in coding regions per acc code
	#
	save_prot_coding_introns($dbh, %prot_id_to_coding_exon_loci, @prot_exon_stats);


	#
	#	Intron / exon counts, intron sizes, disruptions per protein
	#
	save_prot_exon_stats($dbh, @prot_exon_stats);
	close $FH_exonic_mismatch_errors;


	print STDERR "\n";
}


#		
#	Merge disruption data into prot_exon_stats
#
merge_disruptions_to_prot_exon_stats($dbh);



#		
#	saves the number of coding exons per gene
#
save_coding_exon_counts_per_gene($dbh);

$dbh->do("ANALYSE taxon.prot_coding_exons;");
$dbh->do("ANALYSE taxon.prot_coding_introns;");
$dbh->do("ANALYSE taxon.prot_exon_stats");
$dbh->do("ANALYSE taxon.prot_coding_structures");


print STDERR "\tCompleted\n";


log_pipeline_stage($dbh, 'NULL', 23, 'Get sequences and identifiers', $start_time);
$dbh->disconnect();

#remove zero size (empty error/status) files
(-z $_) && unlink $_ for (glob("$dir_pipeline_errors/*"));

